Soumya D. Sanyal

In this post, I'm going to try and see if I can elicit any relationships between the explanatory features and the explained features (expenditure).

This part of the project serves to aid in getting to actually know the data, and in hypothesis making, hypothesis refinement, and model selection; it essentially lays the foundation for the "story" that the project will tell.

I find this stage of the modeling process to be one of the most interesting. It's hard to know if and how your priors about the problem you're looking at will influence the modeling process - I think they almost certainly do. At the same time, there's a great deal of dramatic tension - will my hunch about the problem be supported by the data or will the data show me that it's completely false? Will another insight rise to take it's place? Will there be any story at all? All of these outcomes are entirely possible at this stage, and it's quite nerve wracking.

In any case, I can do my best by examining the data for as much supporting (or refuting) evidence in as objective a manner as possible.

Here are some of the questions I'm hoping to answer at the end of this post:

1. What explanatory features are predictive (if any)?
2. What explains the statistical variation in the dataset? What does the randomness look like? Is there a pattern? Surely it's not all just white noise!
3. Which features are mutually correlated? Do features combine to create stronger or weaker relationships?

Next, we've accumulated a lot of boilerplate code that I need to load to start up my modeling session at this point. Rather than subject you to a lot of scrolling, I'm simply going to save it to a python module called blogging.py, and import it in one step. If you'd like to take a look at the environment I'm importing to build the analysis in this post, it's here.

The next step in my analysis is to explore the data a bit more. I want to see if I can get any more insight into what the underlying data generating process is; the more I can rely on a clean look at the data, the better I'll be able to dispel my own incorrect assumptions about what is or should be happening.

One problem I do have to contend with in working with this dataset is that it comes as weighted survey data. Working with weighted data can be quite tricky. In the very first instance, since the unweighted data does not (by definition) accurately represent the population of interest, it's impossible to draw valid inferences from the standard exploratory (or even regression) procedure without accounting for the weights.

When constructing summary data, as in a histogram, or constructing a linear function of the data, this is not much of a problem. Any linear function $f(x)$ of the data is just as easy to compute as $f(Wx)$, where $W$ is a weight matrix, and variances will just be quadratic in $W$ as expected. However, there are myriad other challenges in working with weighted data; there's a relatively recent discussion of some of these challenges in this paper by Andrew Gelman.

For my immediate purposes - my Insight project - the challenge was to come up with some rough and ready ways to accommodate these weights. While the histogram methods in the Python ecosystem have options for specifying item level weights, not all EDA function calls do. In particular, scatterplots don't (as far as I know). Most model estimation calls that I looked into before starting this project did have some accommodations for item weights, so I figured that if I could get past the EDA stage, I might be fine.

I needed to find a way of getting a sense of the pairwise correlations in the dataset while accounting for the weights. In the end, I found that matplotlib hexbin plots had some support for item weights, and I figured that adding the weighted sum of points to hexbin cells would do well enough without my having to write my own visualization functions. In the end, this is what I constructed:

Let's try this out on a an example pair of features.

In this graph, greater (weighted) mass in any one hexbin results in a higher color count (according to the colorbar scale on the right). More tightly packed data results in a greater plot density.

There looks to be a general tendency for older people to spend more on office based services in this picture, but my sense of the trend is obscured by the outliers above 100k. Of these outliers, most are in the 70-80 year range, with one below 20 years. There's not enough data for me to make a robust case about whether age correlates well with the incidence of really high expenditure, so I'm going to try looking at this again after cutting off the outliers.

This is helpful. It looks like a lot of weighted mass (the aqua region at the bottom) is concentrated near zero expenditure. However, it highlights the fact that expenditure trends higher with age (as does variance in expenditure).

It looks like this is a fairly useful visualization technique as far as quick and dirty goes.

Here are the other groups of correlations I want to inspect:

1. Income
2. Utilization
3. Health conditions
4. Duration without insurance
5. Age of diagnosis
6. BMI
7. Demographic information
8. Expenditures.

I also want to look at the four expenditure categories at the same time. Here's a quick method to help me do that:

It looks like the same pattern holds for the other expenditure categories as for office expenditure when it comes to age.

Looking at income:

Here you get the result that while the level of spending moves away from zero as income rises, the variance of expenditure in all four categories declines with income. An important question is whether the reduced variation in spending at higher incomes is simply a result of having fewer people at those income levels, or if there is a systematic effect (more effective health management, better bargaining power).

Looking at how utilization correlates with expenditure:

I get the following out of these plots:

1. Office based utilization and expenditure are positively correlated. Office based utilization is negatively correlated with all other expenditure.
2. Outpatient utilization doesn't strongly correlate with outpatient expenditure. Outpatient utilization correlates somewhat negatively with other categories of expenditure.
3. ER visits correlate negatively with other expenditures, and not particularly at all with ER expenditure.
4. Inpatient visits correspond positively with inpatient utilization, and negatively with all other expenditure.

Next, I'll look at chronic health diagnoses. There are over 40 in my list, so it'll be useful to break them into groups.

Before we take a look at these plots, I want to quickly mention a very counterintuitive finding that I came across when I looked at these. For many of these chronic conditions, the subpopulation that does not have a positive diagnosis has a much heavier right tail than the subpopulation with the chronic condition.

On the face of it, it looked like people overall were at greater risk of spending more if they did not have a chronic illness than if they did. The paradox disappears when you zoom in on these plots around the first $\$10,000$ spent; it shows that many more people without a chronic condition are clustered at $\$0$ spent, relatively speaking. It's interesting to conjecture what may cause the thinner tails for chronic sufferers - one explanation may just be that the larger subpopulation is naturally accompanied by a larger variance.

The plots below shows the scatterplots to this limited scale for emphasis.

Conclusions:

1. A cancer diagnosis overall drives expenditure for office based services and inpatient services, and has less of a clear effect on outpatient spending and ER spending. For these latter two categories, there is greater variance in spending for the subpopulation without a positive cancer diagnosis.
2. There is very little mass in the graphs measuring spending against the various subdiagnoses of cancer. It seems likely that there will be a corresponding lack of statistical significance of these features in any analysis I do, and it may be well to leave those out and just use the overall cancer diagnosis feature.

Conclusions:

1. Diagnoses of emphysema, bronchitis, and asthma drive spending on office services and inpatient services. They have less of an effect on outpatient and ER services. The right tails in these latter categories is heavier in the negative group.
2. The other features do not seem to be as predictive of increased spending.

Conclusions:

1. Diagnoses of high blood pressure, high cholesterol, coronary heart disease, angina, heart attack, other heart disease, stroke, and diabetes drive spending on office services and inpatient services, while having a limited effect on the other two categories of spending, in which they are associated with thinner right tails.
2. A multiple diagnosis of high blood pressure is probably not very predictive.

Conclusions:

1. Diagnoses of joint pain and arthritis drive spending on office services and inpatient services, while having a limited effect on the other two categories of spending, in which they are associated with thinner right tails.

Conclusions:

1. Diagnoses of ADD/ADHD, limitations in physical functioning, and incidence of pregnancy drive spending on office services and inpatient services, while having a limited effect on the other two categories of spending, in which they are associated with thinner right tails.

Conclusion:

1. It doesn't look like spending is correlated with time elapsed without health insurance.

Conclusion:

1. There doesn't seem to be a strong correlation between the age of diagnosis for any of these chronic conditions and annual spending on any category of healthcare.

Conclusions:

1. No features except insurance status generate significant differences in outpatient expenditure.
2. ER spending is affected by race, marital status, employment status, self-employment status, holding multiple jobs, union status, whether you are offered health insurance by your employer, your overall insurance status (notably whether you have private insurance or a public option), whether your household uses an FSA, where you live (notably, there is more zero clustering for ER spending in the South, and less in the Northeast), and educational attainment.
3. Both office and inpatient spending are affected by race, marital status, employment status, having more than one job, being self employed, union status, whether health insurance is offered through employment, whether public or private options for health insurance are used, whether FSA's are used, insurance status, census region, and educational attainment.

Finally, a look at whether body mass index affects spending:

From looking at this, there doesn't seem to be a strong correlation between spending and BMI in any category. I find this somewhat counterintuitive, since obesity is supposed to be correlated with poor cardiovascular health. I'm inclined to leave this feature in to see if there are any interesting interaction effects that may turn up in the modeling process.

At this point, I've been able to look at correlations between the explanatory features and spending.

Something else I'd like to do before I move on to making modeling decisions is to examine the pairwise correlations between the explanatory features. Apart from giving me a better view of the dataset and patterns within, I find it useful to get a sense of which features provide "orthogonal" signal. This is generally useful in building models that are based on least squares optimization, since parameter estimates tend to exhibit poor precision in the presence of mutually correlated regressors. While this can often be mitigated by regularizing the objective function (for example, penalizing the objective function by a scalar multiple of the $L_1$ norm of the parameter vector will generally send "redundant" parameter estimates to zero), it's nice to know why this might be happening behind the scenes.

A problem to take care of before hand is that almost none of the categorical features are coded in a way where it makes sense to compute a raw correlation coefficient. For example, the levels here:

don't correlate linearly with the severity of the high cholesterol condition.

To more accurately get a sense of the correlations between the properties measured by these categorical features, I'm going to create new dummy variables (that take values either 0 or 1) that correspond to these properties.

I'll do this by groups of predictors:

This should allow me to get a decent sense of the correlations within the explanatory variables (both within these groups and overall).

Here's a heatmap of the correlations between the health status variables.

Conclusions:

1. The arthritis, joint pain, diabetes, limitations in physical functioning and heart health features are mutually correlated. Presumably this arises from comorbidities in advanced age.
2. The limitations in physical functioning feature is correlated with all chronic conditions ranging from cancer diagnosis to high cholesterol, bronchitis, high blood pressure and coronary heart disease.
3. High cholesterol is strongly correlated to the incidence of high blood pressure.

While there are more interesting facts to be gleaned from this heatmap, one thing that strikes me is that the limitations in physical functioning feature is correlated with most chronic conditions. This may be useful to include in a model as a controlling feature in order to reduce omitted variables bias and best tease out the individual contributions to spending of the distinct chronic conditions. As a quick check to see if this might be a problem, let's see if the limitations in physical functioning feature correlates with spending in even one category:

While there's much less mass in the positive bucket, there seems to be much more zero clustering in the negative bucket. For the moment, I might suspect a possibility of introducing omitted variables bias if I were to leave out the limitations in physical functioning feature.

Next, here's how some of the demographic information correlates:

1. There's a natural block of negative correlation between the various race categories - most people self-identify with only one of these, so any self identification with one category reduces the level of self-identification with the others.
2. Being married seems to correlate with education level. Possibly this is just a consequence of greater age.

Let's look at the employment variables:

Conclusions:

1. The occupation sectors most likely to offer health insurance seem to be education/health/social services and manufacturing.
2. People whose employers offer health insurance are most likely to be insured for the whole year.
3. People who have private insurance are likely to have it for the whole year.

Next, a look at what we get if we toss in employment, demographic and health variables into one heatmap:

This plot is quite rich in detail. Some conclusions:

1. Higher educational attainment correlates with employers offering health insurance.
2. High blood pressure, arthritis and heart disease correlate with being on public insurance. Age is likely the intermediating factor.
3. Professional, health/social/education, management/business/financial and manufacturing occupations tend to offer employer insurance.
4. Higher educational attainment correlates with likelihood of being employed.

Finally, a quick look at correlations between diagnosis ages for chronic conditions:

Conclusions:

1. The age of diagnosis for each pair of (high cholesterol, high blood pressure), (high cholesterol, arthritis) and (heart attack, coronoary heart disease) are positively correlated.

To conclude this part, I've learned something about the structure of the explanatory features, and what it will mean for the design matrix in my model. Also, I think heatmaps are a really nice visualization tool.

At this point I have a better view of the data, and a sense of what features I can expect to be predictive. This is a decent foundation for a pretty rich modeling exercise.

In the next post, I'll write about the model selection process.

In this post, I'm going to look at the explanatory features in the dataset.

First, boilerplate:

Let's get our dataset, variables and functions from before.

Next I'm going to prototype a function that will spit out the empirical distributions of each of the explanatory variables.

I'm going to have to treat categorical variables differently to continuous variables.

One rough and ready way to break out the categorical variables is to check how many levels that feature has in the dataset - if there are relatively few (say below a threshold of 10 or 20), I'll assume the feature is categorical, and otherwise I'll assume it's continuous. To deal with some of the edge cases, it's easy to note that anything involving an age, or a # of visits should be continuous.

Let's see if this works.

Let's see if this gave us what we wanted.

Pretty good!

Next I'll prototype a function to produce bar plots for the categorical variables. Since the data points are weighted, I need to make sure I scale each point out correctly before I aggregate under each label as well.

Here's the function:

The categorical features above split up nicely into two groups that are interesting. The first consists of health status features, and the second consists of demographic and socioeconomic features. I'll group them accordingly:

Here are a few observations about these health status features:

1. High blood pressure, diabetes, high cholesterol, joint pain and arthritis are relatively common in the population.
2. Heart disease (angina, heart attack, etc) and stroke are less prevalent, but still substantial.
3. Ditto for lung disease (emphysema, bronchitis, asthma).
4. While a cancer diagnosis is relatively common, the individual cancer diagnoses have very low mass in the "yes" bucket. This will likely pose a problem for estimating models, unless we either upsample, or otherwise penalize the other responses. Including these features in our model may be informative if different cancers affect expenditure in different ways, but we'll have to balance this against the lack of reliability of parameter estimates that we'll get from these features. Unless we have a good reason to believe that different cancers cost very different amounts to treat, we may be best off using the aggregate cancer diagnosis feature, and dropping the individual diagnosis features.
5. Practically no-one uses a peak flow meter.
6. There's relatively incidence of ADD/ADHD in the population.
7. Relatively few people are pregnant at the time of the survey.
8. General limitation in physical functioning is fairly prevalent.

Next let's look at the demographic and socio-economic features.

About the demographic/socio-economic features:

1. The male/female subpopulations are evenly balanced in the population, with slightly more women than men.
2. The race/ethnicity distribution is consistent with other surveys.
3. The distribution estimates that about 70m poeple are under 16, a similar number has never been married, about 130m people are married.
4. About 150m people were employed when the survey was conducted. About 100m people are not employed (not in the labor force, or unemployed). This is consistent with reported data elsewhere. About 20m people are self employed, and about 10m people hold more than one job.
5. Remarkably few people - less than 50m - have a choice of health plans at work. I'm interested in whether people who have a choice of health plans tend to pay less for healthcare on average, assuming there is any effect coming from encouraging people to think about their purchasing choices on subsequent utilization or on prices paid/negotiated.
6. The modal entry for industry group and occupational group (by far) is inapplicable. This is a constructed feature in the dataset, so I assume this wasn't directly asked of respondents. It's hard to say what this might mean, or if it'll be useful to keep in our model.
7. About 10m people are in a union. It would also be interesting to know whether being part of a union affects utilization or expenditure in any way.
8. About 60m people did not hold health insurance at the time. These numbers are consistent with pre-ACA uninsured levels.
9. About 80m people were offered health insurance through their employers, and 40m were not.
10. About 180m people under 65 had private insurance coverage during the survey, about 50m people under 65 had only public insurance options, and about 50m people under 65 were uninsured. There were about 20m people over 65 holding a mix of private insurance and medicare, and a similar number with only medicare coverage.
11. About 30m people lived in a household where someone had a flexible spending account. This is about $10\%$ of the population, but FSAs (especially coupled with high-deductible plans) are an interesting demand-side innovation and it would be interesting to understand if this plays a role in determining expenditure.
12. A negligible number of people (relatively speaking) had coverage under a state insurance plan.
13. The distribution of people over census regions is approximately as follows: about 120m in the south, about 75m people in the west, about 65m people in the midwest, and about 60m people in the northeast.
14. The distibution over educational attainment is as follows: about 70m people have an associate's degree past high school, about 65m have a GED or high school diploma, a little under 60m have an 8th grade qualification or less, about 30m people stopped high school between 9th and 12th grade and do not have a high school diploma, a little over 40m people have a 4-year degree, about 25m people have a masters, doctoral or professional degree, and about 25m people fall into the inapplicable bucket.

Now that we've had a look at the categorical features, let's take a look at the continuous features.

There are broadly speaking, continuous variables in the following categories: ages, income levels, time since some event, number of occurences of an event, and total expenditures.

Let's group these variables:

Let's look at the ages first.

The documentation only mentions that the value $-1$ is recorded for age when there is no valid data collected and no imputation can reasonably be made. In the following graph I'll censor this level, since there's not really a good way to interpret this that I can think of. I will probably also leave it out when building the model.

The first thing to notice is the huge mass at 85. There aren't actually over 6m people in the country of age exactly 85; the data is top-coded at this value to preserve respondent privacy. This means I won't be able to build a model that discriminates on age past 85, and I'll just have to lump everyone older than that into the same bucket. I'll hopefully be able to get a reasonable cross section of that subpopulation using other features.

Next let's look at the age of diagnosis features. Since most people don't have chronic illnesses of these types, the mass at -1 dominates the distribution. Accordingly, I'm going to censor those points in these plots, so each of these distributions are conditional on having that chronic condition.

Here's a modified function from the last post for producing histograms with weighted data.

1. The "lifestyle" diseases: heart attack, angina, high blood pressure, coronary heart disease, stroke, emphysema, high cholesterol, and diabetes all follow very similar distributions, which are all not-quite normal, with modes in the 50-60 range, and with some right skew.
2. It will be interesting to look at correlations between these lifestyle disease features. I don't know a priori if this should be high or low, but I would imagine that having multiple comorbidities would be predictive of higher expenditure on healthcare.
3. Arthritis follows a similar distribution to the lifestyle diseases.
4. Asthma is very heavily skewed left, and is most diagnosed at very young ages, with the diagnosis rate trailing off as the population gets older.
5. The age of diagnosis for ADD/ADHD is very tightly clustered in the 0-20 range, with a bell shaped distribution, non-existent tails, and some skew to the left.
6. A generally interesting hypothesis about the relationship between expenditure and the age of diagnosis of any of these diseases may be that the earlier the diagnosis, the better the patient is at making lifestyle changes that can mitigate expenditure levels. It might be interesting to see if this holds up in the modeling.

Let's look at BMI.

The BMI (body mass index) is computed by dividing a person's weight in pounds by their squared height in square inches and multiplying by 703.

Here's how to interpret the adult body mass index according to the CDC guidelines for interpretation.

The distribution of BMI in the population is estimated by this histogram to be bell-shaped, with a heavier right tail, and centered around 25 or so.

Next, on to the income variable:

What's really weird about this is that I see responses below 0; far below 0!

Let's take a closer look.

Thankfully there are not many data points here. An important question is how to interpret an income of $-258220$. The smaller values for negative incomes are also interesting - are these reported after something like adjusting for public assistance by the respondents?

Let's take a look at the survey documentation.

FAMINC13 contains total family income for each person’s CPS family. Family income was derived by constructing person-level total income comprising annual earnings from wages, salaries, bonuses, tips, commissions; business and farm gains and losses; unemployment and workers’ compensation; interest and dividends; alimony, child support, and other private cash transfers; private pensions, IRA withdrawals, social security, and veterans payments; supplemental security income and cash welfare payments from public assistance, Temporary Assistance for Needy Families, and related programs; gains or losses from estates, trusts, partnerships, S corporations, rent, and royalties; and a small amount of “other” income. Person-level income excluded tax refunds and capital gains. Person-level income totals were then summed over family members, as defined by CPSFAMID, to yield CPS family-level total income (FAMINC13). 

Looking at the expenditure features, there's nothing particularly outlandish about these records apart from the incomes. Given this, I think I'm inclined to go by: "when in doubt, throw it out" when it comes to these data points.

Let's truncate at 0 and take a look at the shape of the distribution.

The distribution is unimodal (mode around $\$50,000$), and with most of the mass concentrated between $\$0$ and $\$200,000$, and a thin (but very well off!) tail past that point.

Let's move on to the durations variables.

Since the majority of respondents did carry insurance in the previous year, I'll drop the points at -1 to get a sense of the distribution conditional on not being insured.

This is quite interesting. Here's what the documentation has to say about this feature:

For persons who were covered by health insurance on January 1st, it was ascertained if they were ever without health insurance in the previous year (NOINSBEF). The number of weeks/months without health insurance was also ascertained (NOINSTM, NOINUNIT). 

I'm interested to know if this feature might be predictive of expenditure. It's hard to say - going a long time without insurance might be reflective of being unemployed, or if the respondent did have the means to buy insurance, may be a proxy for their (lack of) demand for risk mitigation. If either of these correlations hold, this might say something about the effect of either bargaining position, or risk tolerance, on propensity to spend.

Next, the utilization features:

It seems clear that most people use relatively little healthcare. Let's confirm if a lot of these utilizations are at 0.

That's definitely the case. Let's see what the distributions look like conditional on positive utilization.

These distributions all look very similar; the modal positive utilization level is one visit, and utilization rapidly decays, although there are very heavy tails for office based services and inpatient services (we've noticed this already in the expenditure features for these categories).

The main difference between these categories is the degree of use. For example, the number of people in the modal office visit level is over 50m, for outpatient and ER is about 30m, and for inpatient nights is about 5m.

Similar pictures in this case.

We've already looked at the patterns in expenditure, so this takes us to the end of the continuous variables we want to look at.

The next thing I'd like to do is to take a closer look at the features to check for bivariate relationships. Do specific diagnoses determine higher levels of expenditure? Do people who make more money pay more for healthcare? If there is such a relationship, is it linear? Is the growth faster or slower than linear? Do certain diagnoses, while not individually affecting expenditure much, combine to have a large effect, perhaps larger than the sum of their individual pieces?

Looking more closely at the data will, I hope, help me validate or refute interesting hypotheses like these, and possibly generate many more. The first step is to try and visualize the possible relationships, which will be in the next post.

In the last post, I managed to construct the MEPS dataset as a csv file. Here's the result, for reference. (Warning: the file is large: about 250MB.)

The next thing to do is to take a look at the dataset, and see if we need to do any more processing. Since we have it written to disk as a csv file, I can import it using pandas' read_csv function. I may have to do a bit of plotting when I explore the dataset, so I'll go ahead and load all the packages I imagine I'll need for this.

Here's the boilerplate:

Let's load the dataset as a pandas dataframe and see what it contains:

First let's take a look at what features we have:

The first thing I get from this is that I need to learn what these features actually are, so I'll have to look up the MEPS survey methodology or documentation in order to figure out what these features actually are. There are 1792 features for each record, so this may be a bit involved.

The second thing I get out of looking at the column names is that they really do need a bit of cleanup. For example, the first feature is called "DWELLING UNIT ID", but the second is called " PERSON NUMBER" (note the extra space at the beginning). This seems to be a running problem in the column names. There's also a column right at the end with no name.

There are 1792 column names to fix, which is entirely too many to fix by hand.

On the other hand, it looks like the cleanup consists of a couple of simple rules:

Remove the initial space in the feature name
Pop the last entry from the column names

Let's write a function to do the first job:

Let's try this on our dataframe and see if it does the first task.

So far so good; now I just need to get rid of the last column.

Great, that worked. Alright, let's roll that into the clean_columns function:

Perfect.

One thing worth remarking on is that I've actually extracted the variable description rather than the variable name in this step. This was a deliberate decision, as I find it much more intuitive to work with variable names like "MOST RECENT COLONOSCOPY (>39) - R5/3" rather than the variable name "CLNTST53". Those of you familiar with the constraints of variable names in SQL databases will object that this is going to cause trouble in the future. I decided to go with the more descriptive variable names, and write a coding-decoding function to help move data in and out of SQL databases a bit later in the project. Overall, I found this to be a reasonable choice.

To this end, let me go ahead and get rid of the spaces in the variables names before that causes any bugs down the road.

Let's run the data columns through this function:

Great!

Next, I'd like to learn a bit more about what these features are. The main thing to keep in mind is that having a large number of features will make our model very unwieldy.

The crucial thing for this model is for it to be as useful as possible, and if I have to collect 1700 features or so about a person before I can give them any useful feedback, this model won't be used by anyone.

Secondly, for the purposes of my project, I'm most concerned with how the demographic, economic, or health status of a consumer drives cost, since these are the things that can inform how we choose to design the market for healthcare services.

To this end, let's take a look at what's available at the MEPS page that points to the dataset. This page gives us two helpful pieces of information: the first is the documentation for the dataset and the second is the codebook, which we've already used in constructing the csv file.

While the codebook contains the entire listing of variables in the dataset and the start-end points in the ASCII file, it does not tell us what values the variables actually take. To be fair, the variable names are fairly descriptive so we might be able to make educated guesses here, but a quick query of one of the columns in the dataset we've constructed:

raises the question: what do the values of this variable mean? It has a max of 2, and a min of -9. That doesn't correspond to any prior I currently have over the values of a diagnosis of high cholesterol. I'd expect to see things like positive, negative, not tested, and so on. Let's take a more detailed look:

The most frequently populated values are 2, 1 and -1, but there are also records with values -8, -7 and -9.

This warrants a closer look at the documentation. A quick search for "cholesterol" turns up:

High Cholesterol

Questions about high cholesterol were asked of persons aged 18 or older. Consequently, persons aged 17 or younger were coded as “Inapplicable” (-1) on these variables. CHOLDX ascertained whether the person had ever been diagnosed as having high cholesterol. Through 2007, a person-level variable (CHLAGE) indicated the age of diagnosis for high cholesterol on the Person-Level Use PUF. The age of diagnosis for high cholesterol (CHOLAGED) is included in this file. This variable is top-coded to 85 years of age.

which is somewhat helpful, but doesn't explain the use of the values -8 and -9.

Further digging into the documentation yields the following table:

2.2 Reserved Codes

The following reserved code values are used:

Value     Definition
-1 INAPPLICABLE     Question was not asked due to skip pattern
-2 DETERMINED IN PREVIOUS ROUND     Question was not asked in round because there was no change in current main job since previous round
-7 REFUSED     Question was asked and respondent refused to answer question
-8 DK     Question was asked and respondent did not know answer
-9 NOT ASCERTAINED     Interviewer did not record the data
-10 HOURLY WAGE >= $76.96     Hourly wage was top-coded for confidentiality
-13 INITIAL WAGE IMPUTED     Hourly wage was previously imputed so an updated wage is not included in this file

This is helpful; it indicates that the most relevant levels for our analysis ought to be the following three:

1. Positive diagnosis of high cholesterol
2. Negative diagnosis of high cholesterol
3. Inapplicable to respondent

and given the low population in levels below -1, we can either aggregate these data points with level==-1, or discard them entirely. In the following analysis, we choose to discard them.

A final question that needs to be asked at this step is which level corresponds to a positive diagnosis of high cholesterol? A close look at the documentation indicates that level==1 corresponds to having the condition, and level==2 corresponds to not having the condition. Knowing which level is which doesn't particularly affect the computations that generate our models going forward, but it is important for being able to interpret the models.

As a quick sanity check, let's get a visualization of the distribution of the data between these levels:

Of the respondents for whom the question was applicable, it looks like the ratio of level==1 to level==2 is 7/17 or about 2/5. This is close to the CDC estimates of the prevalence of high cholesterol amongst American adults: http://www.cdc.gov/dhdsp/data_statistics/fact_sheets/fs_cholesterol.htm .

Next, instead of using level==1 or level==2 to describe the values that a categorical feature takes, I'd like to translate the levels of these categorical variables into their descriptions so my plots are easier to read.

Here's a dictionary that implements this using the documentation quoted above:

Let's write a function that implements this as necessary for a feature.

The features whose levels I need to interpret in this way are all categorical features; the documentation indicates that the continuous variables aren't generally coded (apart from top-coding to preserve privacy). Accordingly, it would help to be able to automatically determine whether the feature I'm looking at is categorical or continuous. I'll try implementing the following rough rule of thumb: if a feature has less than 10-20 levels, I'll assume it's categorical, and otherwise I'll assume it's continuous.

Let's classify the features in the dataset:

Finally, I'll write a function that interprets the levels of the features in any categorical variable:

Here's an example of how it transforms the levels in a plot:

That visualization was a nice way of learning something about the data, and confirming a hunch we had about how to interpret it. In fact, simple visual plots of different features in the dataset are a crucial way of getting some initial insight into the data, and helping us prioritize the most fruitful approaches for modeling the data.

To fix ideas, let's scope down what we want to go into our model.

Explicitly, we'd like to be able to answer the following questions: which (if any) of the following features about a person determine what they will pay for healthcare in a year?

1. Demographic information:
   1. Age
   2. Sex
   3. Race/Ethnicity
   4. Marital Status
   5. Educational attainment
   6. Census region
2. Income information (these are all top-coded in the dataset for confidentiality reasons):
3. Idiosyncratic health status variables:
   1. High blood pressure, including multiple diagnoses
   2. Heart disease (including coronary heart disease, angina, myocardial infarction, and other unspecified heart disease)
   3. Stroke
   4. Emphysema
   5. Chronic Bronchitis
   6. High cholesterol, including the age of diagnosis
   7. Cancer
   8. Diabetes
   9. Joint pain
   10. Arthritis
   11. Asthma
   12. Attention Deficit Hyperactivity Disorder/Attention Deficit Disorder (ADHD/ADD)
4. Adult Body Mass Index
5. Health Insurance Variables
   1. Public Insurance Indicators
   2. Private Insurance Indicators
   3. Any Insurance Indicators
   4. HMO plan
6. Utilization, expenditure and source of payments information

It's fair to say that this is a lot of information to process. One thing that stands out here is that there are a great number of features already recorded for us in this dataset; this is in sharp contrast to some projects I've done using text and/or social media data, in which the hard part of the modeling process is feature construction in some appropriate way.

Since there is such a rich variety of features, I imagine most of the feature engineering in this project will amount to bucketing variables and levels in appropriate ways to optimize explained variation, rather than constructing new features from scratch.

To this end, let's make a list of the exact features we want to start working with.

1. AGE AS OF 12/31/13 (EDITED/IMPUTED)
2. SEX
3. RACE/ETHNICITY (EDITED/IMPUTED)
4. MARITAL STATUS-12/31/13 (EDITED/IMPUTED)
5. FAMILY'S TOTAL INCOME
6. ADULT BODY MASS INDEX (>17) - RD 5/3
7. HIGH BLOOD PRESSURE DIAG (>17)
8. CORONARY HRT DISEASE DIAG (>17)
9. AGE OF DIAGNOSIS-HIGH BLOOD PRESSURE
10. MULT DIAG HIGH BLOOD PRESS (>17)
11. AGE OF DIAGNOSIS-CORONARY HEART DISEASE
12. ANGINA DIAGNOSIS (>17)
13. AGE OF DIAGNOSIS-ANGINA
14. HEART ATTACK (MI) DIAG (>17)
15. AGE OF DIAGNOSIS-HEART ATTACK(MI)
16. OTHER HEART DISEASE DIAG (>17)
17. AGE OF DIAGNOSIS-OTHER HEART DISEASE
18. STROKE DIAGNOSIS (>17)
19. AGE OF DIAGNOSIS-STROKE
20. EMPHYSEMA DIAGNOSIS (>17)
21. AGE OF DIAGNOSIS-EMPHYSEMA
22. CHRONC BRONCHITS LAST 12 MTHS (>17)-R3/1
23. CHRONC BRONCHITS LAST 12 MTHS (>17)-R5/3
24. HIGH CHOLESTEROL DIAGNOSIS (>17)
25. AGE OF DIAGNOSIS-HIGH CHOLESTEROL
26. CANCER DIAGNOSIS (>17)
27. CANCER DIAGNOSED - BLADDER (>17)
28. CANCER DIAGNOSED - BREAST (>17)
29. CANCER DIAGNOSED - CERVIX (>17)
30. CANCER DIAGNOSED - COLON (>17)
31. CANCER DIAGNOSED - LUNG (>17)
32. CANCER DIAGNOSED - LYMPHOMA (>17)
33. CANCER DIAGNOSED - MELANOMA (>17)
34. CANCER DIAGNOSED - OTHER (>17)
35. CANCER DIAGNOSED - PROSTATE (>17)
36. CANCER DIAGNOSED - SKIN-NONMELANO (>17)
37. CANCER DIAGNOSED-SKIN-UNKNOWN TYPE (>17)
38. CANCER DIAGNOSED - UTERUS (>17)
39. DIABETES DIAGNOSIS (>17)
40. AGE OF DIAGNOSIS-DIABETES
41. JOINT PAIN LAST 12 MONTHS (>17) - RD 3/1
42. JOINT PAIN LAST 12 MONTHS (>17) - RD 5/3
43. ARTHRITIS DIAGNOSIS (>17)
44. AGE OF DIAGNOSIS-ARTHRITIS
45. ASTHMA DIAGNOSIS
46. AGE OF DIAGNOSIS-ASTHMA
47. DOES PERSON STILL HAVE ASTHMA-RD3/1
48. DOES PERSON STILL HAVE ASTHMA - RD 5/3
49. ASTHMA ATTACK LAST 12 MOS - RD3/1
50. USED ACUTE PRES INHALER LAST 3 MOS-RD5/3
51. USED>3ACUTE CN PRES INH LAST 3 MOS-RD5/3
52. EVER USED PREV DAILY ASTHMA MEDS -RD 5/3
53. NOW TAKE PREV DAILY ASTHMA MEDS - RD 5/3
54. HAVE PEAK FLOW METER AT HOME - RD 5/3
55. EVER USED PEAK FLOW METER - RD 5/3
56. ADHDADD DIAGNOSIS (5-17)
57. AGE OF DIAGNOSIS-ADHD/ADD
58. PREGNANT DURING REF PERIOD - RD 3/1
59. LIMITATION IN PHYSICAL FUNCTIONING-RD3/1
60. EMPLOYMENT STATUS RD 3/1
61. HAS MORE THAN ONE JOB RD 3/1 INT DATE
62. SELF-EMPLOYED AT RD 3/1 CMJ
63. CHOICE OF HEALTH PLANS AT RD 3/1 CMJ
64. INDUSTRY GROUP RD 3/1 CMJ
65. UNION STATUS AT RD 3/1 CMJ
66. OCCUPATION GROUP RD 3/1 CMJ
67. HEALTH INSUR HELD FROM RD 3/1 CMJ (ED)
68. HEALTH INSUR OFFERED BY RD 3/1 CMJ (ED)
69. EMPLOYER OFFERS HEALTH INS RD 3/1 CMJ
70. FULL YEAR INSURANCE COVERAGE STATUS 2013
71. ANYONE IN RU HAVE FSA - R3/1
72. # WKS/MON WOUT HLTH INS PRV YR-PN 18 ONL
73. PRIVATE INSURANCE ANY TIME IN R5/R3
74. PUBLIC INS ANY TIME IN R5/R3
75. INSURED ANY TIME IN R3/1
76. ANY TIME COVERAGE BY STATE INS - R3/1
77. # OFFICE-BASED PROVIDER VISITS 13
78. # OUTPATIENT DEPT PROVIDER VISITS 13
79. # EMERGENCY ROOM VISITS 13
80. # NIGHTS IN HOSP FOR DISCHARGES 2013
81. FINAL PERSON WEIGHT 2013
82. TOTAL OFFICE-BASED EXP 13
83. TOTAL OUTPATIENT PROVIDER EXP 13
84. TOT HOSP IP FACILITY + DR EXP 13
85. TOTAL ER FACILITY + DR EXP 13
86. CENSUS REGION AS OF 12/31/13
87. EDUCATION RECODE (EDITED)

It's going to be extremely interesting to see if these have explanatory power for the statistical variation in spending on healthcare!

Speaking of which, what are the variables we're trying to predict?

1. TOTAL OFFICE-BASED EXP 13
2. TOTAL OUTPATIENT PROVIDER EXP 13
3. TOT HOSP IP FACILITY + DR EXP 13
4. TOTAL ER FACILITY + DR EXP 13

The dataset breaks out healthcare expenses into these four location based categories: office, outpatient, inpatient and ER.

Another modeling choice that needs to be made is whether to try and predict total expenses, or predict each category of expense individually. In my Insight project, I opted to try and predict them individually, to see if I could shed some light on whether the service location affected prices charged. There's reason to think this may be the case: ER sticker shock is widely reported in the news, more so than office visits.

Finally, it's time to reflect a bit on what we've done here. We've isolated explanatory variables that we hope and expect should be predictive of how much a person is charged for healthcare services in four different categories. That's great, but there's a very large number of predictors: over 80. Since the ultimate goal is to build a usable interface to the model that people can use to get a sense of their expected outlays on healthcare, it makes sense to refine and reduce this list. But how?

One way to think about this is to consider what information people are likely to be comfortable sharing about themselves that could help to predict expenses. Another is to focus on what the most predictive variables are, and then ask users to populate these variables in order to get a prediction.

In my Insight project, I opted to balance both these concerns by ranking the explanatory variables in terms of statistical variation explained, cutting off at a reasonable threshold (somewhere in between 10 and 20 variables, corresponding to 10-20 questions asked of the user), and also culling variables that would be particularly burdensome to the user to have to report.

The next thing to think about is: now that we have explanatory variables and explained variables selected, how do we best go about building a model?

My next step is to examine the data feature by feature, to get a better sense of what it looks like and to pick up hunches for what may explain the variation. Along the way I'm sure I'll develop some more detailed hypotheses about healthcare spending, and discover that some of what I implicitly assume to be true about the healthcare market isn't actually so.

First, let's look at the explained variables:

This is interesting. It looks like upto half of respondents spend less than $200 for office-based services in a year. The maximum though, is concerning: over a hundred thousand dollars? Let's take a closer look at that.

There are 4 people with office based expenses of over 100k. I'd like to take a look to see if these records are bad data, or if they should be kept in the dataset. My prior is that they should be kept in the dataset (the MEPS survey has been running for over 15 years by people who I can only assume know what they're doing). But it doesn't hurt to take a closer look.

Now, I don't want to look at all the extra features in the dataset beyond those that I've decided I care about, so let's go about trimming the dataframe down to contain just these features.

As a matter of good practice, let's write the current dataset to a file, in case we break something and need to come back to this checkpoint.

I've manually written the explanatory variables to a text file located at "/home/soumya/research/insight/insight_project/modeling/code/dsm_exog_.txt". I'm going to read this into memory:

I'm going to make sure there are no duplicates in my list:

For future reference, let's write this dataset to disk so we don't have to do all this over again.

Let's check this worked: